from numpy import *
from numpy import tensordot as td
import itertools

def main():
  global crystal
  crystal = 'fcc'
  orientation = 'fakerand'
  indices=range(12)
  choices=5

  Mtot = 0.0
  numIter=5000.0

  #result of 5000 random orientations
  # 3.06082265184

  #result of 5000 rand() orientations  

  temp=itertools.combinations(indices,choices) #make a list of all the combinations of indices
  slipIndicesList=list(temp) #convert the combinations object to a list

  slipIndicesArray=array(slipIndicesList) #convert the list to an array... probably can skip  this step
  

  for count in range(int(numIter)):
    phi1, theta, phi2=setFCCOrientation(orientation)
    rmat=createRmatBunge(phi1, theta, phi2)
    
    s0,m0=createSlipDirNorm(rmat)
    s0m0=zeros((3,3,12))
    emat = zeros((5,12))
	
	#Declare the matrix outer(s0,m0) for all the slip systems
    for i in range(12):
      s0m0[:,:,i]=outer(s0[:,i],m0[:,i])
	
	#Declare the matrix e as the slip combinations
    for i in range(12):
      emat[0,i] = s0m0[1,1,i]
      emat[1,i] = s0m0[2,2,i]
      emat[2,i] = s0m0[1,2,i]+s0m0[2,1,i]
      emat[3,i] = s0m0[0,1,i]+s0m0[1,0,i]
      emat[4,i] = s0m0[0,2,i]+s0m0[2,0,i]

    uniqueSlipIndices=array([])  

    #Check prospective values to see if they are unique. If so, retain them. 
    for prosList in slipIndicesArray:

      tempE=transpose(array([emat[:,prosList[0]],emat[:,prosList[1]],emat[:,prosList[2]],emat[:,prosList[3]],emat[:,prosList[4]]]))
      #if linalg.matrix_rank(tempE)==choices:
      if matrixrank(tempE)==choices:      
        uniqueSlipIndices=append(uniqueSlipIndices,prosList)
      #if abs(linalg.det(tempE)) >= 1e-10:		
        #uniqueSlipIndices=append(uniqueSlipIndices,prosList)

    uniqueSlipIndices=uniqueSlipIndices.reshape((size(uniqueSlipIndices)/choices,choices))
  
    #Define things for Taylor Analysis
    eps=1.0
    Mmin=1e7

    strainVec=zeros(5)
    #StrainVec in form (e22,e33,e23,e13,e12) with e11 omitted... it is derived quantity
    strainVec[0],strainVec[1]=eps/2.0,-eps  

    for prosList in uniqueSlipIndices:
      tempET=transpose(array([emat[:,prosList[0]],emat[:,prosList[1]],emat[:,prosList[2]],emat[:,prosList[3]],emat[:,prosList[4]]]))
      #tempGam=td(linalg.inv(tempET),strainVec,(1,0))
      tempGam=linalg.solve(tempET,strainVec)
      if sum(abs(tempGam))/eps <= Mmin:
        Mmin=sum(abs(tempGam))/eps

    Mtot+=Mmin
  
  print 'Mavg: ', Mtot/numIter   
  
  
  
  
def matrixrank(A,tol=1e-8):
    s = linalg.svd(A,compute_uv=0)
    return sum( where( s>tol, 1, 0 ) )  
  
def createRmatBunge(phi1, theta, phi2):
  """Returns 3x3 rotation matrix based on rotation angles in Bunge notation"""
  
  if (not(0.0<=phi1<=2*pi) or not(0.0<=theta<=pi) or not(0<=phi2<=2*pi)):
    print 'Angles are out of bounds in createRmatBunge call'
    return
  
  rmat=zeros((3,3))
  s1,c1,st,ct,s2,c2=sin(phi1), cos(phi1), sin(theta), cos(theta), sin(phi2), cos(phi2)
  rmat[0,0],rmat[0,1],rmat[0,2]= c1*c2 - s1*s2*ct, s1*c2 + c1*s2*ct , s2*st
  rmat[1,0],rmat[1,1],rmat[1,2]= -1.0*c1*s2-s1*c2*ct ,-1.0*s1*s2+c1*c2*ct , c2*st 
  rmat[2,0],rmat[2,1],rmat[2,2]= s1*st, -1.0*c1*st, ct
  
  return rmat

def createSlipDirNorm(rmat):
  """
  This function returns the rotated slip direction and slip normal matrices based on rmat.
  Currently, only has fcc slip systems, but can easily be tailored to do BCC or HCP slip systems.
  """

  if crystal=='fcc':
    s0,m0=zeros((3,12)),zeros((3,12))
    r2,r3=1.0/sqrt(2), 1.0/sqrt(3) #declare values used for normalized vectors
    
    #declare slip systmes for FCC material, s0 is slip direction, m0 is slip normal
    s0[0,0],s0[1,0],s0[2,0],m0[0,0],m0[1,0],m0[2,0]= -r2, r2,  0.0,  r3,  r3, r3, 
    s0[0,1],s0[1,1],s0[2,1],m0[0,1],m0[1,1],m0[2,1]= 0.0, -r2, r2,   r3,  r3, r3
    s0[0,2],s0[1,2],s0[2,2],m0[0,2],m0[1,2],m0[2,2]= r2,  0.0, -r2,  r3,  r3, r3
    s0[0,3],s0[1,3],s0[2,3],m0[0,3],m0[1,3],m0[2,3]= -r2, -r2, 0.0, -r3,  r3, r3
    s0[0,4],s0[1,4],s0[2,4],m0[0,4],m0[1,4],m0[2,4]= r2,  0.0, r2,  -r3,  r3, r3
    s0[0,5],s0[1,5],s0[2,5],m0[0,5],m0[1,5],m0[2,5]= 0.0, r2,  -r2, -r3,  r3, r3
    s0[0,6],s0[1,6],s0[2,6],m0[0,6],m0[1,6],m0[2,6]= r2,  -r2, 0.0, -r3, -r3, r3
    s0[0,7],s0[1,7],s0[2,7],m0[0,7],m0[1,7],m0[2,7]= 0.0, r2,  r2,  -r3, -r3, r3
    s0[0,8],s0[1,8],s0[2,8],m0[0,8],m0[1,8],m0[2,8]= -r2, 0.0, -r2, -r3, -r3, r3
    s0[0,9],s0[1,9],s0[2,9],m0[0,9],m0[1,9],m0[2,9]= r2,  r2,  0.0,  r3, -r3, r3
    s0[0,10],s0[1,10],s0[2,10],m0[0,10],m0[1,10],m0[2,10]= -r2, 0.0, r2,   r3, -r3, r3
    s0[0,11],s0[1,11],s0[2,11],m0[0,11],m0[1,11],m0[2,11]= 0.0, -r2, -r2,  r3, -r3, r3
    
    #return the matrices after they have been rotated by rmat
    return td(rmat, s0, (1,0)), td(rmat, m0, (1,0))
    
  elif symmetry=='bcc':
    print 'Have not implemented bcc slip direction and normal yet in createSlipDirNorm'

    return
  elif symmetry=='hcp':
    print 'Have not implemented hcp slip direction and normal yet in createSlipDirNorm'
    return

  else:
    print 'Invalid symmetry type entered in createSlipDirNorm'
    return

def setFCCOrientation(arg):
  """
  Set orientation for FCC crystal. can take in arguments: '100', '110', '111', or 'rand'. 
  It returns angles in radians in Bunge notation, which corresponds to a z,x',z'' rotation. 
  See the Mathematica notebook find-orientations.nb for calculations of FCC special planes of rotation.
  """
  
  if arg=='100':
    return pi/2.0, pi/2.0, 0.0
  elif arg=='110':
    return 3.0*pi/4.0, pi/2.0, pi
  elif arg=='111':
    return 3.0*pi/4.0, 0.9553166181, 0.0
  elif arg=='rand':
    return random.uniform(0,2*pi),arccos(random.uniform(0.0,1.0)),random.uniform(0,2*pi)
  elif arg=='fakerand':
    return random.uniform(0,2*pi),random.uniform(0.0,pi),random.uniform(0,2*pi)
  else:
    print 'Did not feed in argument 100, 110, 111, or rand into setFCCOrientation(), returning 0,0,0'
    return 0.0,0.0,0.0  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
if __name__ == "__main__":
    main()
